﻿using System;
using System.Windows.Forms;
using System.IO;
namespace ESCoordinateConvert
{
    public enum CoordType
    {
        None,
        Beijing54,
        Xian80,
        CGCS2000,
        WGS84
    }

    public enum ProjectionType
    {
        None,
        UTM,
        MERCATOR,
        GAUSS_KRUGER
    }
    public delegate int ProgressDelegate(double complete, IntPtr message, IntPtr arg);
    public class ProgressClass
    {
        private ProgressBar m_pb;
        public ProgressClass(ProgressBar pb)
        {
            m_pb = pb;
        }

        public int ProgressFun(double complete, IntPtr message, IntPtr arg)
        {
            if (m_pb.InvokeRequired)
            {
                m_pb.BeginInvoke(new ProgressDelegate(ProgressFun),complete,message,arg);
            }
            else
            {
                m_pb.Value = (int)(complete * 100);
            }
            return 1;
        }
    }

    public class GdalTransform
    {
        private static readonly string GCS_Beijing_1954 = "GEOGCS[\"GCS_Beijing_1954\",DATUM[\"D_Beijing_1954\",SPHEROID[\"Krasovsky_1940\",6378245.0,298.3]],PRIMEM[\"Greenwich\",0],UNIT[\"Degree\",0.0174532925199433]]";
        private static readonly string GCS_CGCS_2000 = "GEOGCS[\"GCS_CGCS_2000\",DATUM[\"D_CGCS_2000\",SPHEROID[\"CGCS2000\",6378137.0,298.257222101]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]]";
        private static readonly string GCS_Xian_1980 = "GEOGCS[\"GCS_Xian_1980\",DATUM[\"D_Xian_1980\",SPHEROID[\"Xian_1980\",6378140.0,298.257]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]]";
        private static readonly string GCS_WGS_1984 = "GEOGCS[\"GCS_WGS_1984\",DATUM[\"D_WGS_1984\",SPHEROID[\"WGS_1984\",6378137.0,298.257223563]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]]";

        private static readonly string GAUSS_KRUGER = "PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",LATORIGIN],PARAMETER[\"central_meridian\",CENTERLON],PARAMETER[\"scale_factor\",1],PARAMETER[\"false_easting\",EASTOFFSET],PARAMETER[\"false_northing\",0],UNIT[\"Meter\",1]";
        private static readonly string UTM = "PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"False_Easting\",EASTOFFSET],PARAMETER[\"False_Northing\",0],PARAMETER[\"Central_Meridian\",CENTERLON],PARAMETER[\"Scale_Factor\",0.9996],PARAMETER[\"Latitude_Of_Origin\",LATORIGIN],UNIT[\"Meter\",1]";
        private static readonly string Mercator = "PROJECTION[\"Mercator_1SP\"],PARAMETER[\"False_Easting\",EASTOFFSET],PARAMETER[\"False_Northing\",0],PARAMETER[\"Central_Meridian\",CENTERLON],PARAMETER[\"latitude_of_origin\",LATORIGIN],UNIT[\"Meter\",1]";


        /// <summary>
        /// 获取地理坐标投影wkt
        /// </summary>
        /// <param name="coordType">椭球类型</param>
        /// <returns></returns>
        public static string GetGCSWkt(CoordType coordType)
        {
            string wkt = "";
            switch (coordType)
            {
                case CoordType.Beijing54:
                    wkt = GCS_Beijing_1954;
                    break;
                case CoordType.Xian80:
                    wkt = GCS_Xian_1980;
                    break;
                case CoordType.CGCS2000:
                    wkt = GCS_CGCS_2000;
                    break;
                case CoordType.WGS84:
                    wkt = GCS_WGS_1984;
                    break;
            }
            return wkt;
        }

        /// <summary>
        /// 获取指定坐标系统指定带号的wkt
        /// </summary>
        /// <param name="coordType">椭球类型</param>
        /// <param name="zoneType">分带类型</param>
        /// <param name="zoneNum">带号</param>
        /// <returns></returns>
        public static string GetWktByZone(CoordType coordType, ProjectionType projType, double zoneType, int zoneNum, bool hasZoneNum)
        {
            string wkt = "PROJCS[";
            double centerLon = 0.0;
            string projstr = "";
            switch (projType)
            {
                case ProjectionType.GAUSS_KRUGER:
                    if (zoneType == 3.0)
                        centerLon = zoneType * zoneNum;
                    else if (zoneType == 6.0)
                        centerLon = zoneType * zoneNum - 3;
                    projstr = GAUSS_KRUGER + "]";
                    break;
                case ProjectionType.MERCATOR: //标准MERCATOR东向偏移固定为0
                    projstr = Mercator.Replace("EASTOFFSET", "0") + "]";
                    break;
                case ProjectionType.UTM:  //标准utm东向偏移固定为500000
                    centerLon = (zoneNum - 31) * 6 + 3;
                    projstr = UTM.Replace("EASTOFFSET", "500000") + "]";
                    break;
            }
            projstr = projstr.Replace("CENTERLON", centerLon.ToString()).Replace("LATORIGIN", "0");
            if (hasZoneNum)
            {
                switch (coordType)
                {
                    case CoordType.Beijing54:
                        wkt += "\"Beijing_1954_" + zoneType + "_Degree_GK_Zone_" + zoneNum + "\",";
                        wkt += GCS_Beijing_1954 + "," + projstr.Replace("EASTOFFSET", zoneNum + "500000");
                        break;
                    case CoordType.Xian80:
                        wkt += "\"Xian_1980_" + zoneType + "_Degree_GK_Zone_" + zoneNum + "\",";
                        wkt += GCS_Xian_1980 + "," + projstr.Replace("EASTOFFSET", zoneNum + "500000");
                        break;
                    case CoordType.CGCS2000:
                        wkt += "\"CGCS2000_" + zoneType + "_Degree_GK_Zone_" + zoneNum + "\",";
                        wkt += GCS_CGCS_2000 + "," + projstr.Replace("EASTOFFSET", zoneNum + "500000");
                        break;
                    case CoordType.WGS84:
                        wkt += "\"WGS_1984_UTM_Zone_" + zoneNum + "N\",";
                        wkt += GCS_WGS_1984 + "," + projstr.Replace("EASTOFFSET", zoneNum + "500000");
                        break;
                }
            }
            else
            {
                wkt = GetWktByCustom(coordType, projType, zoneType, centerLon, 0.0, 500000.0);
            }
            return wkt;
        }



        /// <summary>
        /// 获取指定坐标系统自定义投影的wkt
        /// </summary>
        /// <param name="coordType">椭球类型</param>
        /// <param name="zoneType">分带类型</param>       
        /// <param name="centerLon">中央经线</param>
        /// <param name="orginLat">基准纬度</param>
        /// <param name="eastOffset">东向偏移</param>
        /// <returns></returns>
        public static string GetWktByCustom(CoordType coordType, ProjectionType projType, double zoneType, double centerLon, double orginLat, double eastOffset)
        {
            string wkt = "PROJCS[";
            switch (coordType)
            {
                case CoordType.Beijing54:
                    wkt += "\"Beijing_1954_" + zoneType + "_Degree_GK_CM_" + (int)centerLon + "E\",";
                    wkt += GCS_Beijing_1954 + ",";
                    break;
                case CoordType.Xian80:
                    wkt += "\"Xian_1980_" + zoneType + "_Degree_GK_CM_" + (int)centerLon + "E\",";
                    wkt += GCS_Xian_1980 + ",";
                    break;
                case CoordType.CGCS2000:
                    wkt += "\"CGCS2000_" + zoneType + "_Degree_GK_CM_" + (int)centerLon + "E\",";
                    wkt += GCS_CGCS_2000 + ",";
                    break;
                case CoordType.WGS84:
                    wkt += "\"WGS_1984 / " + (int)centerLon + "E\",";
                    wkt += GCS_WGS_1984 + ",";
                    break;
            }
            switch (projType)
            {
                case ProjectionType.GAUSS_KRUGER:
                    wkt += GAUSS_KRUGER.Replace("CENTERLON", centerLon.ToString()).Replace("LATORIGIN", orginLat.ToString()).Replace("EASTOFFSET", eastOffset.ToString()) + "]";
                    break;
                case ProjectionType.MERCATOR:
                    wkt += Mercator.Replace("CENTERLON", centerLon.ToString()).Replace("LATORIGIN", orginLat.ToString()).Replace("EASTOFFSET", eastOffset.ToString()) + "]";
                    break;
                case ProjectionType.UTM:
                    wkt += UTM.Replace("CENTERLON", centerLon.ToString()).Replace("LATORIGIN", orginLat.ToString()).Replace("EASTOFFSET", eastOffset.ToString()) + "]";
                    break;
            }
            return wkt;
        }

        public static bool TransformTiff(string srcfile, string desfile, string wkt, ProgressClass pgc)
        {
            string logpath = Path.GetDirectoryName(srcfile) + "\\err.txt";
            try
            {
                OSGeo.GDAL.Gdal.AllRegister();
                OSGeo.GDAL.Gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");                
                OSGeo.GDAL.Driver drv = OSGeo.GDAL.Gdal.GetDriverByName("GTiff");
                OSGeo.GDAL.Dataset ds = OSGeo.GDAL.Gdal.Open(srcfile, OSGeo.GDAL.Access.GA_ReadOnly);

                //double[] adfGeoTransform = new double[6];
                //ds.GetGeoTransform(adfGeoTransform);
                //adfGeoTransform[0]  左上角x坐标   
                //adfGeoTransform[1]  东西方向分辨率  
                //adfGeoTransform[2]  旋转角度, 0表示图像 "北方朝上"  
                //adfGeoTransform[3]  左上角y坐标   
                //adfGeoTransform[4]  旋转角度, 0表示图像 "北方朝上"  
                //adfGeoTransform[5]  南北方向分辨率

                OSGeo.OSR.SpatialReference src = new OSGeo.OSR.SpatialReference("");
                src.ImportFromWkt(ref wkt);

                OSGeo.OSR.SpatialReference dst = new OSGeo.OSR.SpatialReference("");
                dst.CopyGeogCSFrom(src);
                string dst_wkt;
                dst.ExportToWkt(out dst_wkt);
                src.ExportToWkt(out wkt);

                OSGeo.GDAL.Dataset vrt = OSGeo.GDAL.Gdal.AutoCreateWarpedVRT(ds, wkt, dst_wkt, OSGeo.GDAL.ResampleAlg.GRA_Bilinear, 0);
                OSGeo.GDAL.Dataset dso = drv.CreateCopy(desfile, vrt, 0, null, (pgc != null) ? new OSGeo.GDAL.Gdal.GDALProgressFuncDelegate(pgc.ProgressFun) : null, "transform");

                dso.FlushCache();
                ds.Dispose();
                vrt.Dispose();
                drv.Dispose();
                dso.Dispose();
                src.Dispose();
                dst.Dispose();
                if (File.Exists(srcfile + "esearth.ovr"))
                    File.Move(srcfile + "esearth.ovr", srcfile + ".ovr");
            }
            catch(Exception ee)
            {
                File.AppendAllText(logpath, ee.Message);
                return false;
            }
            return true;
        }

        public static bool TransformAsc(string srcfile, string desfile, string wkt, ProgressClass pgc)
        {
            string logpath = Path.GetDirectoryName(srcfile) + "\\err.txt";
            try
            {                
                OSGeo.GDAL.Gdal.AllRegister();
                OSGeo.GDAL.Gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
                OSGeo.GDAL.Driver drv = OSGeo.GDAL.Gdal.GetDriverByName("AAIGrid");
                OSGeo.GDAL.Dataset ds = OSGeo.GDAL.Gdal.Open(srcfile, OSGeo.GDAL.Access.GA_ReadOnly);
                OSGeo.OSR.SpatialReference src = new OSGeo.OSR.SpatialReference("");
                src.ImportFromWkt(ref wkt);
                OSGeo.OSR.SpatialReference dst = new OSGeo.OSR.SpatialReference("");
                dst.CopyGeogCSFrom(src);
                string dst_wkt;
                dst.ExportToWkt(out dst_wkt);
                src.ExportToWkt(out wkt);

                OSGeo.GDAL.Dataset vrt = OSGeo.GDAL.Gdal.AutoCreateWarpedVRT(ds, wkt, dst_wkt, OSGeo.GDAL.ResampleAlg.GRA_NearestNeighbour, 0);
                OSGeo.GDAL.Dataset dso = drv.CreateCopy(desfile, vrt, 1, null, (pgc != null) ? new OSGeo.GDAL.Gdal.GDALProgressFuncDelegate(pgc.ProgressFun) : null, "transformasc");
             
                dso.FlushCache();
                ds.Dispose();
                vrt.Dispose();
                drv.Dispose();
                dso.Dispose();
                src.Dispose();
                dst.Dispose();
                if (File.Exists(srcfile + "esearth.ovr"))
                    File.Move(srcfile + "esearth.ovr", srcfile + ".ovr");
            }
            catch (Exception ee)
            {
                File.AppendAllText(logpath, ee.Message);
                return false;
            }
            return true;
        }

        public static double GetLonLatResolution(string srcfile, string wkt)
        {
            try
            {
                OSGeo.GDAL.Gdal.AllRegister();
                OSGeo.GDAL.Gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
                OSGeo.GDAL.Dataset ds = OSGeo.GDAL.Gdal.Open(srcfile, OSGeo.GDAL.Access.GA_ReadOnly);

                OSGeo.OSR.SpatialReference src = new OSGeo.OSR.SpatialReference("");
                src.ImportFromWkt(ref wkt);
                OSGeo.OSR.SpatialReference dst = new OSGeo.OSR.SpatialReference("");
                dst.CopyGeogCSFrom(src);
                string dst_wkt;
                dst.ExportToWkt(out dst_wkt);
                src.ExportToWkt(out wkt);

                OSGeo.GDAL.Dataset vrt = OSGeo.GDAL.Gdal.AutoCreateWarpedVRT(ds, wkt, dst_wkt, OSGeo.GDAL.ResampleAlg.GRA_NearestNeighbour, 0);
                double[] adfGeoTransform = new double[6];
                vrt.GetGeoTransform(adfGeoTransform);

                //adfGeoTransform[0]  左上角x坐标   
                //adfGeoTransform[1]  东西方向分辨率  
                //adfGeoTransform[2]  旋转角度, 0表示图像 "北方朝上"  
                //adfGeoTransform[3]  左上角y坐标   
                //adfGeoTransform[4]  旋转角度, 0表示图像 "北方朝上"  
                //adfGeoTransform[5]  南北方向分辨率

                ds.Dispose();
                vrt.Dispose();
                src.Dispose();
                dst.Dispose();
                return adfGeoTransform[1];
            }
            catch
            {
                return 0.0;
            }
        }
    }
}

